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Abstract 

In this article we investigate the feasibihty of constructing stable, local bases for com- 
puting with kernels. In particular, we are interested in constructing families b = (fo^)^gH 
that function as bases for kernel spaces S{k,E) = {X1{gh ^) I (o?)4eH G so that 

each basis function can be obtained by very few kernels 

6^ = ^ A^^^k{-,^) Ai;,^ = for all but a few ^. 

This is reminiscent of the construction of the B-spline basis from the family of truncated 
power functions. 

We demonstrate that for a large class of kernels (the Sobolev kernels as well as many 
kernels of poly harmonic and related type) such bases exist . In fact, the basis elements 
can be constructed using a combination of roughly C'(logiV)'^ kernels, where d is the local 
dimension of the manifold and N is the dimension of the kernel space (i.e. N = #S). 
Viewing this as a preprocessing step - the construction of the basis has computational 
cost O [N (log N)"^) . Furthermore, we prove that the new basis is Lp stable and satisfies 
polynomial decay estimates that are stationary with respect to the density of S. 

1 Introduction 

The purpose of this article is to investigate robust bases for spaces associated with a positive 
definite or conditionally positive definite kernel /c : M x M — )• M, where M is a C°° (closed) 
compact Riemannian manifold. The dimension of M is d. The kernels that we discuss below 
belong to a wide class that includes the thin-plate splines and similar kernels when M = 
or SO{3). 
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The spaces associated with a kernel from this class are defined as follows. Let H C M be 
a finite set of points, called centers, having cardinality N = ^E. The centers are scattered in 
the sense that they do not need to belong to a regular grid. In the positive definite case, the 
space S{k, E) associated with the kernel k and the set E is just S{k, E) = span{A;(-, E), ^ G E}; 
that is. 



S{k,S) 



The conditionally positive definite case is similar; we will discuss it it in Section [4.11 - specif- 
ically in (j4.3p . For these spaces, if h := max^g= dist(x, ^) is the mesh norm (fill distance), 
q := ^ min^^^ dist(T/, ^) is the separation radius and ps '■= h/q is the mesh ratio, then as long 
as />= < POi where po is fixed, Lebesgue constants are uniformly bounded and approximation 
rates for functions in Sobolev spaces W^{Wl) are 0{h"^), with the constants independent of 
other properties of E [201 ttH] • 

Two other remarkable properties of S{k, E) concern its Lagrange basis, {xc(")}cgs- Recall 
that in a Lagrange basis each basis function satisfies Xii'n) = ^i,'q when r] E E. What was 
shown in [20^ [19] is that decays exponentially fast away from ^ for special kernels, and 
algebraically fast for many others. 

Equally as important, as we shall prove below in Theorem 14.31 if we express the x^'s in 
the standard basis, 

where the coefficients A^^^ are well-known to be the entries of the inverse of the interpolation 
matrix, then |j4g,r;| decays as a function of dist(7/,.^) at the same rate as |xc(3^)l decays in 
dist(x, ^) - i.e., exponentially or algebraically, as the case may be. Prior to our work, the only 
provable results concerning decay of these coefficients were done by Fornberg [T3] in the case 
of M and for gridded data, using Fourier techniques that do not carry over to the scattered 
case. 

The difficulty with the Lagrange basis is that each X£, is computationally costly both to 
construct (as a linear combination of A;(-,^), ^ G H) and to compute with. Are there better 
bases? Here is what we would desire in a basis for S{k, E). 

Each basis function should be highly localized and nearly scalable with respect to the 
mesh norm h of E. By this we mean that each basis element is of the form 

where the ??'s come from small subset of the centers E and satisfy the following requirements: 

i) / 0} = c(#H) 

/ dist(x, 



bp(x)\ < a 



h 



where the cost c(A^) is constant or slowly growing with N = ^E and the function cr(x) 
decays rapidly: at an exponential rate cr(x) < Ce"*^'^' or at least at a fast polynomial rate 
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a{x) < C(l + The B-spline basis, constucted from the family of truncated power 

functions, is a model kernel method and in the d = 1 case provides an ideal solution to the 
problem we consider. 

The main results of this paper demonstrate that such kernel bases exist and that each 
basis function can be computed in nearly fixed time. Viz., an individual basis function can 
be computed in ©((log A^)'^) time while the full basis can be computed in 0{N{log N)'^) time 
with = #H = Ch~'^. Moreover the basis is stable. 

The main tool employed is Theorem 14.31 which allows one to bound the rate of decay of 
Lagrange coefficients in terms of the corresponding decay rate of the Lagrange functions. In 
particular if the Lagrange functions have exponential decay so too do their corresponding 
coefficients. While, as we mentioned earlier, this fact was previously known in the scaled 
lattice case, but no such estimates have been available in the scattered case. 

The sphere As an example of our main results for the sphere and the restricted 
surface splines of order s + 1, given by ks+i{x, a) := (l — x-a)^ log(l — x-a), for s = 1, 2, 3, ... , 
we have the following theorem, which is a corollary of Theorem 15.11 in Section [5l 

Theorem. For a sufficiently dense set of centers and for a sufficiently large constant r 
there is basis (b^)^^s for the space S'(A;, H) satisfying the following. 

• Each basis element b^ = A^^(^k{-, Q) is composed of at most M := r(log A^)^ kernels. 
I.e. 

c{#E) = M = T{logNf . 

• Each basis element exhibits polynomial decay: there exist constant C and J for which 

Mx)\<c{i+^^y\ 

• The rate of polynomial decay J depends linearly on the constant of proportionality r by 

j = oiV^). 



• The basis is Lp stable: there are constants < ci < C2 < oo depending on t so that for 
all sequences a = (0^)5^= £ I^" the following holds: 



cig2/P||a||^^(H) < 



Lp{S2) 



The restricted surface splines on are of special importance, largely because the sphere 
is the setting of many problems of scientific interest, but also because the kernels themselves 
have a convenient, closed form representation and their approximation power that is well 
understood and optimal in the sense that approximation rates are in line with smoothness 
assumptions for the target functions (i.e., approximands in Sobolev classes Wp or Besov classes 
-Bp 00 are approximated by functions in S{k,E) with error decaying like h'^). 
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We give some numerical examples for §^ with such kernels (and others) in Section [6j 
Precondioners Over the years practical implementation of kernel approximation has pro- 
gressed despite the ill-conditioning of kernel bases. This has happened with the help of clever 
numerical techniques like multipole methods O El E] and often with the help of precondi- 
tioners O [121 ISSl [35] of which [71 [271 [2H] sue of special interest to us, because these involve 
attempts to construct local bases. Indeed, another offshoot of our results in that one can now 
estimate, a priori, the number of coefficients needed to guarantee good preconditioners. Many 
results already exist in the RBF literature concerning preconditioners and "better" bases. For 
a good list of references and further discussion, see [11]. Several of these papers use "local 
Lagrange" functions in their efforts to efficiently construct interpolants. The number of points 
chosen to localize the Lagrange functions are ad hoc and seem to be based on experimental 
evidence. For example, Faul and Powell, in [13], devise an algorithm which converges to a 
given RBF interpolant that is based on local Lagrange interpolants using about thirty nearby 
centers. Beatson-Cherrie-Mouat, in [3], use fifty local centers (p. 260, Table 1) in their con- 
struction along with a few "far away" points to control the growth of the local interpolant at 
a distance from the center. In other work. Ling and Kansa [23] and co-workers have studied 
approximate cardinal basis functions based on solving least squares problems. Thus one goal 
of this paper is to provide some theoretical groundwork that may yield future improvements 
in preconditioner algorithms and better bases for kernel spaces. 

Organization We devote Section [2] to treating some pertinent results and definitions 
for Riemannian manifolds. In Section [3] we consider the stable, local bases constructed in 
[201 [TBI [T9] , which have many desirable properties but are computationally infeasible due to 
their cumbersome construction — each basis function of this type requires #H nonzero kernel 
coefficients in its construction and, moreover, computing these requires operations. 
An analysis of these coefficients show that they drop off rapidly — this is demonstrated in 
Section m The rapid decay of these coefficients leads to the (theoretical) existence of efficiently 
constructed bases, but sadly does not indicate the desired construction - this is treated in 
Section [5l In Section [6l we give numerical evidence to bolster the results of the previous 
sections, by giving results of experiments that show how rapidly the Lagrange basis and the 
coefficients decay. In this section we give some examples of techniques that fail to deliver, 
and provide some examples of families that seem to have the desired properties which have 
not been validated theoretically. 

2 Geometric background 

Throughout this paper, M denotes a compact, complete d-dimensional Riemannian manifold. 
The Riemannian metric for M is g, which defines an inner product gp{-, •) = (•, ■)g^p on each 
tangent space TpM; the corresponding norm is | • \g^p. 

The Riemannian metric is employed to measure arc length of a curve 7 via Ijlg^pdt. 
Geodesies are curves 7 : M — )• M that locally minimize the arc length functional giving rise to 
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a distance function ^ 

dist{p,q) = min / Ijlg^pdt. 
7(0)=p Jo 
7{1)=<? 

We denote the ball in M centered at x having radius r by B{x,r). Given a finite set H C M, 
we define its mesh norm (or fill distance) h and the separation radius q to be: 

/i. := sup dist(x, H) and q ■= - inf dist(^,C)- (2.1) 

The mesh norm measures the density of H in M, the separation radius determines the spacing 
of H. The mesh ratio p := h/q measures the uniformity of the distribution of H in M. We 
say that the point set H is quasi-uniformly distributed, or simply that H is quasi-uniform if H 
belongs to a class of finite subsets with mesh ratio bounded by a constant /jq. 

The metric g also induces an invariant volume measure d/i on M. The local form of the 
measure is dfi{x) = ^^/det{g)dx^ ■ ■ ■ dx"^, where det{g) = det{gij). We indicate the measure of 
subsets C M by vol(r2). The integral, and the Lp spaces for 1 < p < oo, are defined with 
respect to this measure. The embeddings 

C(M) C Lp(M) for 1 < p < oo and Lp(M) C -Lg(M) for 1 < g < p < oo 

hold. In addition, L2 is a Hilbert space equipped with the inner product (•,•)• {f ■: d) ^ {f ■• 9) 



Sobolev spaces on subsets of M Sobolev spaces on subsets of a Riemannian manifold 
can be defined in an invariant way, using the covariant derivative (or connection) V (cf. [l]) 
which maps tensor fields of rank j to tensor fields of rank j + 1. The kth. covariant derivative 
of a function is a rank k tensor field and is denoted V'^/. For k = 1, the covariant derivative in 
local coordinates is simply the usual expression for the "gradient" - it can be written simply 
as (V/(a;))j = ^~f{x). For k = 2, the "Hessian" tensor involves Christoffel symbols F^ and 

can be expressed as (V^/(x))ij = -^-[^{x) - I]m=i<^ ^i^C^)^^!^)- Higher order covariant 
derivatives have an analogous expression, using higher order derivatives of the Christoffel 
symbols - see [201 Eqn. (3)]. 

Definition 2.1 ([H p. 32]). Let Q cM. be a measurable subset. We define the Sobolev space 
W^{n) to be all f : M —?■ R such that, for < k < m, |V'^/|g,p in ^2(^2) with associated 
norm 



m,n '■— \\J Ww^in) 



coming from the Sobolev inner product 

{f,9)m,n '■= {f-,9)w^{vt) 



\ 1/2 

j^\V^f\lpd„{p)\ , (2.2) 
:=E / (V'/,^^^ MV). (2.3) 



fc=0 

When = M, we may suppress the domain: {f,g)m = {f,g)m,M dnd 
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Metric equivalence The exponential map allows us to compare the Sobolev norms we've 
just introduced, to standard Euclidean Sobolev norms as follows: 

Lemma 2.2 ([20, Lemma 3.2]). Form S N and < r < im/S, there are constants < ci < C2 
so that for any measurable Q C Br, for all j £ N, j < m, and for any pq S M, the equivalence 

cillnoExpp^ < Mwi{E.p^^^m ^ C2||uoExpp^ H^^.^^^ 

holds for all u : Exppp(il) —5- M. The constants c\ and C2 depend on r and m but they are 
independent of and po . 



3 The Lagrange basis 

For a manifold M, a positive definite kernel A; : M x M — t- M and a set of centers H C M, we are 
concerned with the robustness of the Lagrange basis (xc)cgh for S{k,E), where Xc(0 = ^^,( 
for all ^ S H. The Lagrange basis plays a central role in most interpolation problems, and 
certainly this is the case for radial basis function and kernel interpolation. Decay of the 
Lagrange basis and analytic consequences have notably been considered in [2 ^ 19} [29 } [13]. 



3.1 The kernels considered 

More recently, [201 IISl H] , develop a theory for fast decay and stability of the Lagrange basis 
associated with certain positive definite and conditionally positive definite kernels. 

• "Sobolev kernels" denoted by Km were introduced in [20] for any smooth, complete, 
compact and connected Riemannian manifold. These are the reproducing kernels for 
the Sobolev inner product0for Ty^"(M) when m > d/2: 

m „ 

• This was extended in [19] to treat a broader class of kernels on certain manifolds called 
kernels of polyharmonic and related type (these are discussed in Section [4. 2p . 

• Included in the class considered in [19] are restricted surface splines on S'^ which are 
kernels of the form km{x, a) = (j){x ■ a) where 




(1 - t)™-«f/2 for d odd 

(1 - t)"'-d/2 iog(i _ t) for d even. 



These kernels are conditionally positive definite, meaning that interpolants are con- 
structed by adding an auxiliary function. (In this case, a low degree spherical harmonic.) 
See Section [4T] below. 

^In fact, the inner product can be weighted as X^JILo Cfc /q (V*/, V'^g)^ ^ d/i(p) with non-negative weights 
Ck for which co and Cm are postive - each such reweighting gives a different inner product and a different 
Sobolev spline. 
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The expansion of the functions (p in terms of Gegenbauer polynomials by Baxter and 
Hubbert, [2], leads to Fourier (spherical harmonic) expansions of the kernels, and from 
there to their characterization as Green's functions for elliptic differential operators. 
These operators are of polyharmonic type - they are of the form Q(A) = nj=i(^ ~ ^i) 
for some real numbers ri, . . . ,rm- This, in turn, permits an understanding of the ap- 
proximation power of the kernel, as investigate in [261 I17j : for functions having Lp 
smoothness s up to order 2m (namely, for target functions in smoothness spaces includ- 
ing B^^g(S'^),VFp^(S'^),C7"(S^) with s < 2m), 

distp(/,5(A;,H)) =0(/i^). 

Here, the space S{k^ H) is modified by addition of low degree spherical harmonic terms 
n = {Yi^m I ^ < [m — d/2\} (this is described in Section [^31 below) . 

• Surface splines on 50(3) which are of the form k{x,a) = (j){u:{a~^x)) with ui{x) the 
angle of rotation of x (which is a left and right invariant metric on the group) and 

0(t) = (sin(t/2))'"-^/^ 

In [21], an expansion of (p in even Chebyshev polynomials of the second kind leads to 
a Fourier (Wigner D-function) expansion of the kernel k. As in the spherical case, 
this leads to its characterization as a Green's function for an operator of polyhar- 
monic type on SO{'i), and to a realization of its approximation power: again, for 
/ G BI,^q{SO{2,)),W^{SO{2,)),C'{SO{2,)) with s < 2m we have 

distp(/,5(A:,H)) =0(/i^). 

Restricted kernels An alternative approach, taken in [15], is to consider the manifold 
M as embedded in an ambient Euclidean space , and to use the restriction of a radial basis 
function - a Euclidean (conditionally) positive definite kernel satisfying rotational symmetry 
(of which there are many prominent examples) - as a (conditionally) positive definite kernel 
on M. In a sense, this is a completely different approach, in the sense that such kernels are 
almost never fundamental solutions to differential operators, a key point of jl9j . On the other 
hand, such kernels may be easily localized in the ambient space R", which may lead to an 
effective way of localizing and preconditioning the restricted kernels. Although the theory 
developed in Sections 3 and 4 does not address such kernels, we include a numerical example 
in Section [6l 

3.2 Analytic properties of the Lagrange basis 

The theory developed in [MlllEllin] addresses analytic properties of bases for S{k,E), related 
to locality, stability of approximation and interpolation. In particular the following are shown. 
Locality. The Lagrange basis is a local bases for S{k,'B). That is, 

teW|<Cexp(-.^i!iMl). 
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Stability of interpolation. Interpolation is stable: the Lebesgue constant is bounded 
(and more generally, is the p norm of the interpolant is controlled by the ip norm of the data) . 
Lp conditioning. There are constants depending only on ci,C2 such that ci||a||^p < 

II SjLi ogX? Il-Lp ^ C2||a||£p, with ci,C2 depending only on m, M and the mesh ratio p. In 
particular, they are independent of #H = dim(5'(fc, H)), and, after a suitable normalization, 
independent of p. 

Marcinkiewicz-Zygmund property. The space S{k, H) possess a Marcinkiewicz-Zygmund 
property relating samples to the size of the function. For s G H), this means that the 
norms ||^ i— )• s(0||^ (h) and \\s\\l^ are equivalent, with constants involved independent of 

Stability of approximation in Lp. Approximation by L2 projection is stable in Lp for 
1 < p < 00. In particular, the orthogonal projector with range S'(/c,H) can be continuously 
extended to each Lp, and it has bounded operator norm independent of #H. 

4 Lagrange function coefficients 

In this section we give theoretical results for the coefficients in the kernel expansion of Lagrange 
functions. In the first part we give a formula, relating these coefficients to native space inner 
products of the Lagrange functions themselves (this is Proposition 14. 2p . We then obtain 
estimates on the decay of these coefficients for a class of kernels on certain compact Riemannian 
manifolds (two point homogeneous spaces). 

4.1 Interpolation with conditionally positive definite kernels 

The kernels we consider in this article are conditionally positive definite on the compact 
Riemannian manifold. As a reference on this topic, we suggest |10i Section 4]. 

Definition 4.1. A kernel is conditionally positive definite with respect to a finite dimensional 
space n if, for any set of centers the matrix Ch := {k{C, C))^ ^-^ positive definite on the 
subspace of all vectors a G satisfying ^^^^"^^(O = for p G 11. 

This is a very general definition which we will make concrete in the next subsections. 
Given a complete orthonormal basis {(pj)jen, of continuous functions (i.e., ||(/'j||oo = 1) any 
kernel 

kix,y) ■.= ^k{j)ipj{x)ipj{y) 

with coefficients k G ^2(N) for which all but finitely many coefficients k{j) are positive (neg- 
ative) is conditionally positive definite with respect to IIj- = span(0j | j E J), where 
J = {j \ k{j) < 0}, since, evidently, 
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provided a^<pj{(,) = for j satisfying k{j) < 0. 

In this case if the set of centers H C M is unisolvent with respect to Hj- = span(93j | j G J') 
(meaning that p S Hj and p{^) = for ^ S H implies that p = 0) then the system of equations 

has a unique solution in x for each data sequence y = G C^. 

By writing the same system in matrix form, with collocation matrix Kh = {k{^: C)) ^ ^g=2 
and auxiliary matrix <I> = j)g=x:7' 

:) (:) ^ 

When data is sampled from a continuous function at points S (i.e., = /(C)) that are 
unisolvenlH for lij this solution generates a continuous interpolant: 

/b/ = 4,:7,h/ = J^a5A;(-,0+P/ 

where = X^jej ^jVi ^ and o-^PiC) = for all p € Hj. Indeed, this interpolant is 

unique among functions from the space 



S{k, E) := S{k, E, Uj) := \ a^k{;0 +Pf 



PfeUj, Y,^iPiO = o\ (4.3) 



It has a dual role as the minimizer of the semi- norm ||| • ||| j- induced in the usual way from the 
"native space" semi-inner product 



When u,v £ S{k,E) - meaning that they have the expansion u = ai,g^(-, C) + Pu 

and V = a2,^^('; + Pv with coefficients (aj^^)gg= _L (11^)1= for j = 1,2 - then the 

semi-inner product is 

{u, v)k^j = ^Y1 «i,5«27A;(C, C) 

We can use this expression of the inner product to investigate the kernel expansion of the 
Lagrange function. 



^Meaning that p|h = for p £ Hj implies that p = 0. 
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Proposition 4.2. Let k = X^jeN ^{3)'-Pj'Wj ^ conditionally positive definite kernel with 
respect to the space Hj = spauj^j ipj , and let H be unisolvent for Hj. Then Xri £ S{k,'E) 
(the Lagrange function centered at rj) has the kernel expansion Xrjix) = X^^gs +PC 

with coefficients 

Proof. Select two centers C,r] £ E with corresponding Lagrange functions X( and Xrj £ 'S'(^) S)- 
Because and are both orthogonal to (nj^)!^, we have 

Now define P := <I)(cI)*(I))~i<I)* : £2(S) (11^)1= C ^2(S) to be the orthogonal projection 
onto the subspace of samples of Hj- on H and let P-^ = Id — P be its complement. Then for 
any data y, ()4.ip yields coefficient vectors A and b satisfying P-'"A = A and P-^^h = 0, 
hence P-^K^P'^A = P-'-K^A = P-^y. Because P^ : £2(S) — ^ ^2(S) is also an orthogonal 
projector, and therefore self-adjoint, it follows that 

{xd^),Xvix))k,j = (KhA^,A^)^2(s) = (KhA^,P^A^)^2(h) = {P^KsA(;,Arj)e,(E) 
= (P^ef,A^)^2(H). 

In the last line, we have introduced the sequence = ((5(^,^)5gH for which K^A^ +Pcls = 
which implies that P-'-K^A^ = P-'-e^^. Using once more the fact that P^ is self-adjoint, and 
that A^ is in its range, we have 

{xd^),Xri(.x))k,j = {P^e^,An) = (e^,P-^A^) = (e^, A^) 

and the lemma follows. □ 

4.2 Estimating Lagrange function coefficients 

In 120^119] . it has been shown that Lagrange functions decay rapidly away from the center. We 
can use this characterization of the Lagrange function to estimate the decay of its coefficients. 
In this section we use Proposition 14.21 to estimate the size of coefficients first for the class of 
strictly positive definite functions developed in j20j . Then we attempt to do the same for the 
more general class of kernels of polyharmonic and related type of |19| . 

Sobolev kernels on compact Riemannian manifolds We begin by considering kernels 
k = Km with native space inner product given by an expression like 

•= <^^''^)«™,M / I3{u,v)xdfi (4.5) 
Jm 

where /3 is a pointwise bilinear form fi{u,v)x '■= Y^=Q^j^''^^^^'^)x with the conditioiJl that 
m > d/2 and co,Cm / and cj > for all j = 0, . . . ,m. Such kernels were considered in 

^in such cases, the kernel Km is naturally (strictly) positive definite and moreover, Km is the fundamental 
solution to the elliptic operator £m ~ "^^^o 
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|20j and existence was demonstrated for all d-dimensional, connected, compact Riemannian 
manifolds. 

In this case, the inner product (|4.5p is a Sobolev inner product, and it has a natural gen- 
eralization to inner products for subsets namely {u,v)^^ ^ = {u,v)^rmi^Q'j = l3{u,v)xdiJ,. 
It is possible to estimate KxCjXc)^™!- 

Kx5,xc)fcml < Wxd ^nn (M\b(g, -^'^^t^' o 1 1 xcl I Wo"' (M) + llxdlw2™(M)llxcllyi/m(M\b(c,^i24^)) 

By way of |2U' Corollary 4.4] we have that 

Unfortunately, this family of kernels is not suitable for treating practical problems. In 
particular, the kernels having native space inner products of the form ()4.5p . even when M is 
the sphere, are not known to have closed form representations in terms of the spatial variable 
(despite being zonal and having a simple and satisfying Fourier-Legendre expansion). 

To remedy this, we remove the restriction that the coefficients Cj are non-negative (al- 
though Cm must be positive). An immediate consequence of this is that we must contend with 
a conditionally positive definite kernel. The upshot is that, for a large class of interesting 
manifolds (including spheres and projective spaces) we can write the Dirichlet form (|4.5p as 
linear combinations of powers of the Laplace-Beltrami operators. The motivation for this 
approach is that the restricted surface splines are fundamental solutions for operators of this 
type. We now describe this. 

Polyharmonic and related kernels on 2 point homogeneous spaces Let M be a 

compact, two point homogeneous space. Included among these are spheres, SO{3) and various 
projective spaces. For our purposes, this is a metric space with distance function dist(x,y) 
and measure /i for which /x(b(x,r)) = ^{y \ dist(x,7/) < r} ~ r'^. Because it is compact, there 
is a Laplace-Beltrami operator A with countable spectrum cr(A) = {Aq, Ai, . . . }. Denote the 
corresponding orthogonal, L2 normalized eigenfunctions for A by ('i/'j)jeN- 

For such a manifold and for any A; G N, the operator (V*'')*V'^' can be expressed as 
^^^q6i,A-' with bk = (—1)'^- Consequently, any operator of the form X]j=o 
expressed as Yl^=o^j^'^ with bk = (— l)'^Cfc and vice-versa: 

m m 

V(6o, ■■■hm) 3(co, . . . , c„) with bm = {-iTcm and ^ 6,A^' = ^ c,{V^YVK (4.6) 

j=0 i=l 

Suppose that the kernel km : M x M — t- M acts as the Green's function for the elliptic 
operator Cm '■= Sjlo ^i^"' ~ Q{^)i ™ the sense that 

/=/ km{-,a)Cm[f{a) - pf{a)\(la + pf 
Jm 

where Pf{x) is the projection on 11^;-, pf = YIU L2(M)'^i-, and the complementary part of 
the spectrum of Cm, {Q(Aj) | j ^ J} C a{Cm), is real and lies to one side of (without loss. 
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we can take cr{Cm) C (0, oo) - namely, by considering —k if needed; this is equivalent to taking 
bm > since the spectrum of A has oo as an accumulation point and (Afc)™" > | Z^j^^ ^ji^kVl 
for all but finitely many k). Such a kernel is said to be of polyharmonic or related type. 



The native space "inner product" on subsets It follows directly that km is condition- 
ally positive definite with respect to Hj. What's more, when £^11^ = {0}, the native space 
semi-inner product can be expressed as 



Jm 



with /3{u, v)^ = E^^o Cfc(V^w, V'^v):, and Co ) • • • 1 Cm guaranteed by (j4.6p . The latter expression 
allows us to extend naturally the native space inner product to measurable subsets Q of M. 
Namely, 



This has the desirable property of set additivity: for sets A and B with fi{A (1 B) = 0, we 
have {u, v)AuB,km,J = (^) '^)A,km,J + (^) ^) B,km,J- Unfortunately, since some of the coefficients 
Cfc may be negative, j3{u,u) and {u,u)Q^k,n,j ™ay assume negative values for some u: in other 
words, the bilinear form {u,v) i— )• {u,v)Q^k,n,J is only an indefinite inner product. 

However, when Q has Lipschitz boundary and u has many zeros, we can relate the 
quadratic form Itilif^ fc^ = {u,u)Q^k„i,j to a Sobolev norm UtiH^mj-f^^. Arguing as in [191 
(4.2)], we see that 



If uIh = on a set H with < ho with ho determined only by the boundary of 

(specifically the radius and aperture of an interior cone condition satisfied by d^l), Theorem 
A. 11 of [H] guarantees that < Ch'^\u\w^{Q) with C depending only on the order 

m, the global geometry of M and the roughness of the boundary (in this case, depending only 
on the aperture of the interior cone condition). Thus, by choosing h sufficiently small, h < h* , 
where h* satisfies the two conditions 

h* < ho and C{h*f x (max|c,|) < (4.7) 

j<m 2 



we have 



I l|2 ^ III |||2 ^ j I II II l|2 



The threshold value h* depends on the coefficients cj as well as the radius Rq and aperture 
4>Q of the cone condition for fi. When 0, is an annulus of sufficiently small inner radius, the 
cone parameters can be replaced by a single global constant, and can be taken to depend 
only on Co, ... ,Cm- In other words, only on km ~ cf. [19^ Corollary A. 16]. 
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A direct consequence of this is positive definiteness for such functions, |||tt|||f^ j > with 
equahty only if n|n = 0. From this, we have a version of the Cauchy-Schwarz inequality: if u 
and V share a set of zeros Z (i.e., u\z = v\z = {0}) that is sufficiently dense in O, then 

follows (sufficient density means that h{Z,n) < h* as above). 

Decay of coefficients for kernels of polyharmonic and related type Fortunately, 
Lagrange functions have many zeros, and [19^ Lemma 5.1] guarantees that the Lagrange 
function satisfies the bulk chasing estimate there is a fixed constant < e < 1 so that for 
radii r less than a constant rM depending on M (for a compact, 2-point homogeneous space, the 
injectivity radius is rM = diam(M)/2) the estimate \\Xii\\w^{M\h{i,r)) < ^\\Xi\\w^{M\h{^,r-J2_)) 
holds. In other words, a fraction of 1 — e of the bulk of the tail ||xj ||vK2'"(M\b(5,r)) is to be found 
in the annulus b(^,r) \ b(^,r — ^) of width ^ oc h (with constant of proportionality ^ 
depending only on M, m and the boundary of il). Provided r < im, it is possible to iterate 
this n times for < r. It follows that there is u = — 4/iologe > so that 

IIX?llw2"'(M\b(|,r)) < ^"IIX^IIvi/a^CM) < Ce '^^ ^ \\x^\\w^ (M) ■ 

By |19l (5.1)] (a simple comparison of to a smooth "bump" (p^ of radius q - also an 
interpolant to the delta data {6^), but worse in the sense that lllxglifc j — lll<?^lllfc„ j ~ see the 
proof of Theorem 14.31 below) we have 

\\xdwr{M\hi^,r)) < Cq^l^-^e-^l. (4.9) 
This leads us to our main result. 

Theorem 4.3. Let M 6e a compact, 2-point homogeneous manifold and let km be a kernel 
of polyharmonic or related type, so that the associated elliptic operator Cm annihilates the 
polynomial space lij- Let p > be a fixed mesh ratio. 

There exist constants h* , v and C depending only on M and km if C M. is suf- 
ficiently dense (i.e., /i(H,M) < h*) then the coefficients of the Lagrange function xc, = 
EceH^a^ml-,^ +V(,^ S{E,J) satisfy 

\Au\ < Cq''-^'^ exp (^-zv^MjM) ^ . (4.10) 

Proof. By Proposition 14.21 and set additivity, we have that 

I^CCI = {X£,,xdk^,J = {Xi^X()n(,k^,J + (xc,Xc)n^,fcm,J> 

where 0^ = {a G M | dist(a,C) < dist(a,,^)} , Jig = {a e M | dist(a,,^) < dist(a,C)} , and 
(modulo a set of measure zero) M\ il^ = ilg. For a compact, 2-point homogeneous space, 
and ilg are two balls of radius diam(M)/2. 
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We can apply the Cauchy-Schwarz type inequality ()4.8p to obtain 

Since C b^(C, r) := M \ b (C, ^dist(e, C)) and C h^{C, r) := M \ b idist(e, C)) , we 
can again employ set additivity and positive definiteness (this time |||xcIIIq^ k j — llx^Mk 
which follows from the fact that M = il^^ U 17^ and that vanishes to high order in ilg - the 
same holds for xc) obtain 



(llxcllH^r(b=(c,r)) \lxdkr^,j + \lxdk^,j llx?llH/r{b^(c,r)) 



The full energy of the Lagrange function can be bounded by comparing it to the energy 
of a bump function - for this is which can be defined on the tangent space by using 
a fixed, smooth, radial cutoff function a: (j)^ o Exp^(x) = a{\x\/q). This is done in [191 (5.1)] 
and we have that Ix^ilMfc j ^"^^ i/<^ciMA: j ^'^^ bounded by Cq'^/'^~'^. 

On the other hand, we can employ (14. 9p to treat HxcllWa^C^^CC^)) llx$llw2'"(b''(CA0)' 
which gives 

From this, the result follows. □ 

Note 1. A similar argument shows that, on a compact, 2-point homogeneous manifold, the 
Lagrange function coefficients A^^^ for a general kernel km of polyharmonic and related type 
(regardless of whether Cm annihilates 11 j-) decay like 

1^5,(1 < Cg'^"^™ max(exp(-iy^), h'^'^). (4.11) 

In such cases, the theory developed here and in [19j indicate a slower decay for Lagrange 
functions and coefficients (although it remains an open problem to determine if these rates 
can be improved, and by how much). In particular, this holds for the restricted surface splines 
km{x, a) = {1 — X ■ a)™"'^/^ on odd dimensional spheres (d S 2N + 1 and m > d/2) as well as 
the surface splines on SO (3) (see Section [3]). 

Note 2. Theorem 14 . 3 1 holds for restricted surface splines km{x, a) = (1— x-a)™'"'^/^ log(l— x-a) 
on spheres of even dimension {d € 2N and m > d/2) - in particular for S^. See Section [6] 
below for some numerical examples in this setting. 

5 A better basis: truncating the Lagrange basis 

We now want to show the existence of a good approximation to the Lagrange function xs, 
that uses many fewer elements in its kernel expansion than the N needed for xs, itself. To do 
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this, we win start with the expansion = Y2(s ^5,C^('' approximate it by a truncated 
expansion of the form 

where T(^) C H n B{(^,r{h)). Our goal is to show that, under the assumptions listed below, 
which apply to a wide class of kernels, we may take r{h) = Kh\ log(/i)|, where > ^, while 
maintaining \\x^ — XiWoo < C/i"^, J := Kv — 2m. 

How many basis elements are used in expanding x^? Doing a simple volume estimate 
shows that the number required is 

#T(C) = Oi{Kh\logh\)''/q'') = 0{\loghf) = 0((logiV)^) « N, 

where we have used h/q = p and N = 0{h~''^). 

One final remark before proceeding with the analysis. Finding xs, requires knowing the 
expansion for x^ and carrying out the truncation above. This is expensive, although it does 
have utility in terms of speeding up evaluations for interpolation when the same set of centers 
is to be used repeatedly. The main point is that we now know roughly how many basis 
elements are required to to obtain a good approximation to x^- We are currently engaged 
in investigating cost effective algorithms to obtain approximate Lagrange functions similar to 

First assumptions We make the following three assumptions 

1. The Lagrange functions decay at a rate < C^exp ^—ui^^^^^^^ 

2. The kernel coefficients of the Lagrange function decay like l^^.c;] < CcCxp ^—i'^ ^^^^^'^^ 

3. The Lagrange basis is Lp stable in the sense that 



<C2(?'^/P||a||,^(H). 



We note that the family of restricted surface splines on S'' when d G 2N satisfy these three 
conditions (conditions 1 and 3 are in [19], while condition 2 follows from Theorem 14. 3p . as do 
the Sobolev splines on any compact Riemannian manifold M (condition 1 follows from [20j . 
condition 3 from [18] and condition 2 from Theorem 14.31 again) . 



Decay By the estimate of coefficients in Theorem 14.31 it suffices to retain only the part of H 
that is within Kh\ log h\ from ^, since the coefficients we cut out have size roughly Ch'^~'^^h^^ . 
There are no more than #H < Cd,ph~'^ of them on the d-sphere, and the kernel is uniformly 
bounded, so we have that 

\xi{^)-xd^)\<Ch^''-^^, 
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and we should choose K > ^ least. Indeed, the pointwise estimate above shows that 

which indicates that we may wish to choose K even larger. This is at our discretion, but to 
preserve stability, we choose K > ^n^H^. 

Computational efficiency Since we retain only the coefficients centered at a distance of 
0{h\ log h\) from ^, we use 

#T({) = o('<Mi^Uo((logiV)^) 

coefficients (when centers are quasiuniform) to compute each basis function x^. 

Stability By the L^o stability of the Lagrange basis, for s € S{k,'B), the samples s|s =: 
(A^)^,=s are bounded in the ioo norm by Mis^. Using the same coefficients but in the new 
basis x^, we form s = ^?X$- The difference between the original and new function is 



\s — s 



oo 



< PIkoc E IXe(^) - Xd^)\ < C\\A\U^h^-^-^^q-'' < Ch^'-^-^^-^lU (5.1) 



so 



\s 



loo 



> (1 - > ci (1 - Ch^''-'^"'-'^) (5.2) 

ll^lloo < (l + C/l^"-2"^-'^)||s||oo<C2(l + C/l^''-2"^-'^)P||oo. (5.3) 

This can be viewed in two ways. 

• Provided h is small enough, the family is a basis. There are elements and they are 
linearly independent. (In particular, if the function is zero, all the coefficients are zero.) 

• The family (x^) is stable in L^o, since the map (A^) s is boundedly invertible. 

We stress that it remains to be determined how actually to compute the basis — we have 
simply shown that a preconditioner exists that has complexity O ((log A'^)'^). 

Theorem 5.1. Let Wl be a compact, 2-point homogeneous manifold and let k„i be a kernel 
of polyharmonic or related type, so thai the associated elliptic operator Cm annihilates the 
polynomial space lij- Let p > be a fixed mesh ratio. 

For sufficiently dense S, with h = h{E, M) < H, with H a constant depending only on p, 
M and k„i, there is a basis (&^)^gs whereach basis element = XIi^gs ^C,C^m(") C) composed 
of kernels centered in the ball B{^, —Klogh). The following are satisfied: 

• The cost of constructing each is #{C | 7^ 0} < r (log #2)*^ with r < CK'^. 



16 



• Each basis element exhibits polynomial decay of degree J := Kiy — 2m: there exists C 
for which 



Mx)\ <c{i + 



dist(a;, ^) 



h 



-J 



The basis is Lp stable: there are ci,C2 for which 



<C2<z''/P||a||,^(H). 



ip(M) 



Proof. It remains to demonstrate the Lp stability of (6^) = (x^) for 1 < p < oo. 

When p = 1, we consider a sequence a = (a^)ggH £ ^-iid set s := X^^^Xg- Holder's 
inequality gives 

lis - 5||i,(M) < C||a||,,(=)Vol(M)/i^'^-2- (5.4) 

since for each x we have the estimate \s{x) — s{x)\ < ^X^^^s la^l^ max^^s \x^{^) ~ X^ix)\- 
Interpolating between ()5.4p and ()5.ip (i.e., interpolating the finite rank operator a i— >■ (s — s)) 
gives 

Ho Sll ^ r'?,'f^'^-2m~d(l-l/p)||„|| 

\\s - s\\l^(m) < Cti '-^'||a||^p(H). 

Therefore, 

Noll /^J,-'^''^— 2m~ii(l— 1/p) II _ II ^ ||~|| ^ ||„|| , /^-LKu—2m—d{l — l/p) \\ „\\ 

PllLpiM) - l|a|Up(H) S ||S||Lp(M) S ||S||Li(M) + L-Z^ l|a|Up(H) 

and we have 

ci9"/n|a||,^(=)(l - < ||s||^^(M) < C2g'^/n|a||,^(H)(l + C/i^^-^— '^). 

□ 

Note 3. For a general kernel km of polyharmonic and related type (where £^11^- ^ {0}), the 
estimate for the decay of Lagrange function coefficients A^^i^ is too slowly to guarantee stability 
of the truncated "basis". In this case, we can guarantee only that tail of the coefficients is 
uniformly bounded, \A^^(^\ < Ch'^, and the best estimate we can give to a truncated Lagrange 
function is 1x^(2;) - X^x)] < C. 



6 Numerical examples 

In this section we give some numerical illustrations of the previous results. In the first example, 
we provide results for restricted surface splines on that support Theorem 14.31 In particular, 
we demonstrate that the constants C and which govern the rate of decay, are in fact quite 
reasonable. In the second example, we illustrate how the results from Section [5] can be used 
for practical computations of surface spline interpolants on §^ that involve large point sets. 
Finally, in the last example, we investigate the decay rate of the Lagrange coefficients for the 
restricted surface spline to the Torus, a manifold not covered by the present theory. 
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Geodesic Distance From Center 

Figure 1: Maximum latitudinal values of the Lagrange function for the kernel k2{x,a) = 
{1 — X ■ a) log(l — X ■ a). This experiment was carried out in double precision arithmetic and 
the plateau at roughly 10~^^ occurs due to ill conditioning of the collocation matrices and 
truncation error. 

Example 1. We demonstrate the decay of Lagrange functions and their coefficients for the 
second order restricted surface spline (also known as the thin plate spline) k2{x,a) = (1 — 
X ■ a)log(l — X ■ a). The interpolant takes the form = S^es ^?,C^('' ~^ P^^ where 
is a degree 1 spherical harmonic. In this example, we use the "minimal energy points" of 
Womersley for the sphere - these are described and distributed at the website |36]. The 
value of these point sets is as benchmarks. Each set of centers has a nearly identical mesh 
ratio. Furthermore, the important geometric properties (e.g., fill distance and separation 
distance) are explicitly documented. Their potential theoretic properties and their importance 
in constructing quadrature rules and spherical designs, which are discussed in |321 133j . are 
not pertinent to this work. Because of the nice geometric properties of the minimal energy 
point sets, it is sufficient to consider the Lagrange function centered at the north pole 

e = (0,0,1). 

Figure [U displays the maximal latitudinal valueqj of log^glx^l- We clearly observe the 

*The function is evaluated on a set of points {cj>, 9) with no equispaced latitudes (j) G [0, vr] and ni 
equispaced longitudes 6 £ [0, 2%] 
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Figure 2: Plot of coefficients for a Lagrange function in the kernel space S{k2, 
iment was carried out in double precision arithmetic. 



.This exper- 



exponential decay of the Lagrange functiorH 



\x^{x)\ < C^exp -UL 



guaranteed by [191 Theorem 5.3]. From this figure, the value of ul, which measures the rate 
of exponential decay is observed to be close to 1.35. 

We can visualize the decay of the corresponding coefficients in the same way. We again 
take the Lagrange function centered at the north pole: for each € S, the coefficient of 
the kernel A:(-, (^') in the expansion = ^(,(k{-,C)~^P^ is plotted with horizontal coordinate 
sin((^'). The results for sets of centers of size = 900, 2500 and 10000 are given in Figure 2. 
The exponential decay 



-2m 



exp 



dist(g,C) 
h 



guaranteed by Theorem 14.31 is clearly in force, and we can estimate the constants Uc and 
Cc for the decay of the coefficients following the method used for the Lagrange functions 
themselves - we note that coefficients are shifted vertically, which is a consequence of the 



factor of q' 



d—2m 



q in the estimate (14.10p . Table [T] gives more results with some added 



detail, including estimates of the constants Cl and Cc- 



^At least until a terminal value of roughly 10 at which point there is a plateau beyond which the values 
no longer decay - see below and Figure [3] for an explanation of this. 
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N 


hx 


PX 




Cl 




Cc 


400 


0.1136 


1.2930 


1.1119 


0.8382 


1.0997 


69.9891 


900 


0.0874 


1.5302 


1.3556 


1.0982 


1.3445 


231.5573 


1600 


0.0656 


1.5333 


1.3513 


1.2170 


1.3216 


324.8534 


2500 


0.0522 


1.5278 


1.3345 


0.9618 


1.3117 


470.6483 


5041 


0.0365 


1.5304 


1.3395 


1.1080 


1.3158 


1087.8 


10000 


0.0260 


1.5421 


1.3645 


1.1934 


1.3369 


2564.9 



Table 1: Estimates of fill distance /i, mesh ratio p for some minimum energy point sets on the 
sphere and u and C values for the kernel A;2(x, a) = {1 — x ■ a) log(l — x ■ a). 



The perceived plateau present in the Lagrange function values as well as the coefficients 
shown in Figures [1] and [2] is due purely to round-off error related to the conditioning of 
kernel collocation and evaluation matrices. These results were produced using double-precision 
(approximately 16 digits) floating point arithmetic. To illustrate this point, we plot the decay 
rate of the Lagrange coefficients for the 900 and 1600 point node sets as computed using 
high-precision (40 digits) floating point arithmetic in Figure O The figure clearly shows the 
exponential decay does not plateau and continues as the theory predicts. 

Example 2. In this example we construct a basis (xg)5eH for the kernel space S{k2, H) C C(S^) 
by using ©((logA^)^) centers to construct each Xi- 

With this basis, we use an equivalent representation in the form 

I^f = Y.^M-)^ (6-1) 

where each is a local Lagrange function about the node ^ formed by M ^ basis elements 
of S(A;2,H). Specifically, let T{i) C H such that ^ G T(^), #T(^) = M, and Uges^(0 = 
then 

4 

^? = 2Z H,C.K-^C) + ^h,3'Pr (6-2) 
C6T(0 j=i 

The coefficients a^^^ and 6gj are determined from the conditions 



as: 



The linear system for determining the interpolation coefficients in (jO.ip can be written 

where (KH)i,j = k2{S,i,£,j) and Bij = (t)j{ii), i,j = 1,...,N. The matrix is a A^-by-A^ 
sparse matrix where each column contains M entries corresponding to the values of a^^t^ in 
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Lagrange Coefficient Decay (40 digit aritlimetic) 
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Figure 3: Plot of coefficients for a Lagrange function in the kernel space S{k2,'^)- This 
experiment was carried out in Maple with 40 digit arithmetic. 

(j6.2p . The matrix is a 4-by-A^ matrix with each column containing the values of b^j in 
()6.2p . With the linear system written in this way, one can view the matrix [Ay By]'^ as a 
right preconditioner for the standard kernel interpolation matrix. 

If the sets T(^) are chosen appropriately then the linear system (|6.3p should be "numeri- 
cally nice" in the sense that the matrix KhtIt + ^By should have decaying elements from its 
diagonal and should be well conditioned. In the example below, each T(^) is chosen as the 
M — 1 nearest nodes to ^. Section [5] suggests taking M = 0((logA^)^). Through trial and 
error we found that choosing M = 7[(log;^Q N)"^)] gave very good results over several decades 
of A^. Each set T(^) can be determined O(logA^) operations by using a KD-tree algorithm 
for sorting and searching through the nodes H. The cost for constructing the KD-tree is 
0(A^(log A^)^). Thus, constructing all the sets T takes 0(A^(log A^)^) operations. 

To solve this linear system we will use the generalized minimum residual method (GM- 
RES) [31]. This is a Krylov subspace method which is applicable to non-symmetric linear 
systems and only requires computing matrix-vector products. Ideally, there should be a 
method for computing these matrix vector products in 0{N) or 0{N log N) operations to 
make GMRES more efficient. Keiner et. al. have shown that this can be done in the case of 
the kernel matrix using fast algorithms for spherical Fourier transforms [22]. In the re- 
sults that follow, we have not used this algorithm, but have instead just computed the matrix 
vector products directly. We will investigate the use of these fast algorithms in a follow up 
study. 
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For the numerical tests we use icosahedral node sets H C of increasing cardinality. These 
were chosen because of their popularity in atmospheric fluid dynamics (see, for example, \T6\ 
[3^ [30l I25| ) where interpolation between node sets is often required. The values of / were 
chosen to take on random values from a uniform distribution between [—1, 1]. Table [2] displays 
the number of GMRES iterations to compute an approximate solution to the resulting linear 
systems ()6.3p . As we can see, the number of iterations stays relatively constant as N increases 
and does not appear to increase with A^. 







Number GMRES iterations 


N 


m 


tol = 10-^ 


tol = 10~s 


2562 


84 


7 


5 


10242 


119 


5 


7 


23042 


140 


6 


7 


40962 


154 


5 


7 


92162 


175 


6 


8 


163842 


196 


5 


7 



Table 2: Number of GMRES iterations required for computing an approximate solution to 
()6.3p using icosahedral node sets of cardinality A'^. m corresponds to the number of nodes 
used to construct the local basis and tol refers to the tolerance on the relative residual in 
the GMRES method. The right hand side was set to random values uniformly distributed 
between [—1,1] and the initial guess for GMRES was set equal to the function values. 

Example 3. A second example shows similar results for Lagrange functions for the kernel 
k{x, a) = {x — a)^ log(2; — a) restricted to a torus of outer radius 4 and inner radius 2. In 
other words, the surface parametrized by 

X = (3 + cos v) cos u 
y = (3 + cos v) sin u 
z = sinv, 

with u,v G [0, 27r]. This combination of kernel and manifold is not treated in |19] (the curved 
torus is not even a symmetric space) although it is considered in [TSj , as a subset of that 
is Hi unisolvent. Indeed the torus is the zero set of a degree 4 polynomial in M'^, a sufficiently 
dense subset will also be Hi unisolvent. Hence, interpolation on such sets by shifts of k is well 
posed. 

For this experiment, we consider "minimum energy" point sets H produced by Ayla Gafni, 
Doug Hardin and Ed Saff which have previously been used in [15]. In each case we fix 
the point ^ = (4, 0, 0) and we examine coefficients (^g^^)^gH of the Lagrange function = 
S^GS ^5,C^("' + To solve the interpolation problem, we add a linear polynomial, p £ 
ni(R^), and require the coefficients to satisfy the side conditions J2(eE ^^,C9(.0 = 0, Vg G 
ni(]R^) - in other words, the coefficients annihilate linear polynomials. 

Despite the fact that Theorem 14.31 does not apply, a certain exponential decay is observed 
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for these coefficients as well. As before, the rate of decay seems to be independent of h as 
well as iV = #H. 

We provide two ways to visualize the decay of the coefficients: first by arranging them 
latitudinally with horizontal axis representing the distance from ^ in the u direction and then 
longitudinally with horizontal axis representing the distance from ^ in the u direction. For 
the specific choice of ^ = (4,0,0), both correspond to geodesic distances. (Other choices of ^ 
are observed to have the same rate of coefficient decay.) 




Torus Coefficients 
long trip 



O N=50D 
> N=10D0 
□ N=40D0 




latitude distance from center 



longitude distance from a 



Figure 4: Log plot of coefficients (in absolute value) of Lagrange functions centered at (4, 0, 0) 
for the kernel k{x, a) = \x — a\'^ log |a; — a| using minimal energy sets on the torus of cardinality 
on 500, 1000 and 4000. 

A more complete account of this experiment is given in Table O where estimates of the 
exponential decay rates u^at and VLong (in the v and u directions, respectively) are given. We 
note that they are between 1 and 1.3 in either direction. 



N 


hx 


PX 


^Lat 


Chat 


^Long 


CLong 


500 


0.3383 


1.5242 


1.2312 


6.9946 


1.1504 


7.3231 


750 


0.2737 


1.5024 


1.1556 


9.4737 


1.2122 


15.4474 


1000 


0.2375 


1.5014 


1.2870 


15.6376 


1.2401 


20.4220 


1999 


0.1639 


1.4725 


1.1213 


25.3917 


1.2719 


49.0776 


3000 


0.1333 


1.4479 


1.0793 


36.1593 


1.2421 


58.7687 


4000 


0.1151 


1.4498 


1.1836 


57.4460 


1.2738 


105.2720 



Table 3: Results for a Lagrange function experiment on the torus using the thin plate spline 
kernel k{x, a) = \x — log |x — a|. In addition to fill distance h, mesh ratio p for minimum 
energy set of Saff and Hardin, estimates of the latitudinal and longitudinal exponential decay 
rate and constant z/ and C values are given. 
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